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Executive summary 

Studies of earthquake catalogues for naturally occurring seismicity led to the formulation of the 
Gutenberg-Richter law [1] which States that the probability P(m) that an observed event in a 
sequence of earthquakes has a moment magnitude greater than m satisfies the relation P{m) = 
exp {-b\r\(\()')(m — m 0 )j, whete m 0 is the observation threshold. Evidence has emerged since the 
formulation of the law that the same relation applies for induced seismicity caused by human 
activity, an observation, together with the power-law character of the law when translated into 
moments rather than moment magnitudes, that the Earth’s crust may be in a self-organized 
critical state [2], The use of a robust method for estimating the value of b from seismic data is 
crucial for quantifying the risk of rare but powerful events. In the case of the Groningen field, 
seismic activity due to gas production started in 1991 and has continued ever since. Moment 
magnitude data is provided by KNMI (the Dutch Meteorological Institute) and is rounded to the 
nearest 0.1. This means that the data has to be regarded as being binned, with the bin boundaries 
located at odd multiples of 0.05. 


In this report the maximum likelihood estimate is applied to this binned data and a maximum 
likelihood estimator b MEE is obtained for the whole catalogue, and various temporal and spatial 
subsets. We remark that only events with magnitude greater or equal to 1.5 and occurring after 
the end of April 1995 are considered in this analysis, as only for such events can reliable detection 
be assured [3], The fact that the data bins are of equal width in moment-magnitude space means 
that b MLE can be determined analytically. Monte-Carlo simulations of synthetic catalogues are used 
to determine estimates of the error in b MLE . 


The main conclusions of this study are the following: 

1. The value of b MEE for the entire Groningen catalogue (l st January 1995 — 31 st December 
2014) is consistent with the commonly accepted value b- 1. 

2. While a subdivision of the catalogue into four subsets of equal numbers of events 
suggests more variability in the maximum likelihood estimator for the subsets than would 
be expected if the Avalue of the underlying process were constant, a closer examination, 
involving further subdivisions of the catalogue, and the use of 95% confidence intervals 
obtained using simulated data, does not show any systematic dependence of b value on 
event rate. 

3. Focusing on regions of 5 km radius centered on Loppersum and Ten Boer, we find that, 
while the b MEE for Ten Boer is consistent with b- 1, the b MEE for Loppersum is 
significandy lower, to the extent that we can be 99.9% confident that the b value of the 
underlying process is less than 1. 

4. The temporal sequence of events for each of the regions around Loppersum and Ten 
Boer is split into two equal sub-catalogues and the value of b MEE determined for each sub- 
catalogue. For both the Loppersum and Ten Boer regions, the two values of b MEE were 
found to be consistent with an underlying b that is unchanging in time. 
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1. Introduction 


The moment magnitude m is a logarithmic scale unit used to describe seismic moment M and has 
no physical significance in itself. Other examples of logarithmic scale units are decibels for 
acoustic energy, stellar magnitude for star brightness and gauge for wire diameter. The moment 
magnitude m corresponding to a seismic moment M is given by [4] 

777 = {log l0 (M/M c )}/d (1.1) 

In Eq. (1.1), dis a dimensionless constant with the value 1.5 andAf c = 1.259 x 10 9 N-m. We note 
that Eq. (1.1) can also be written as 

777 = {ln(M / M c )}/{d ln(10)} (1.2) 


According to the Gutenberg-Richter law, the probability P{m) that an observed event in a 
sequence of earthquakes has a moment magnitude greater than m satisfies the relation 

Pijn) = exp{-&ln(10)(7?7 - 777 0 )} (1.3) 

where m 0 is the observation threshold. If M 0 is the threshold seismic moment corresponding to 
we have, from Eq. (1.2), 

77 ? 0 = {ln(M 0 / M c )}/{d ln(10)} (1.4) 

Subtracting Eq. (1.4) from Eq. (1.2), and rearranging the right-hand side yields 

777 - 777 0 = {ln(M / M 0 )}/{d ln(10)} (1.5) 


Substituting Eq. (1.5) into Eq. (1.3) yields the result 

F(M) = P(m) = (M /M 0 )~P (1.6) 

where [3 = b! d (1.7) 


and F(M) denotes the probability that the seismic moment of an observed event is greater than 
M. For induced seismicity, a commonly accepted value of b is unity, corresponding to a [3 value 
of 2/3. The distribution F(M) is a long-tailed Pareto distribution and has infinite variance for < 
2 and infinite mean for /?< 1. Since the threshold for reliable detection of a seismic event can be 
specified, F(M) can be regarded as a one-parameter distribution with parameter (3. We also note 
that Eq. (1.6) allows for M to be unbounded, which is unphysical. For the Groningen field, upper 
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bounds on the maximum possible value of M have been specifïed [3], but such an event is so rare 
that the maximum moment of an observed event in the Groningen catalogue to date is far below 
it. Hence the existence of a “corner moment”, at which F(M) starts to reduce at a rate faster than 
the long tail in Eq. (1.6), cannot be inferred from the data. The focus of this report is the 
estimation of the parameter /? from the observed data (which directly gives b via Eq. 1.7). This is 
done both for the entire catalogue, and also for temporal and spatial subsets. 


Chapter 2 describes the methodology used, which employs the maximum likelihood estimator, 
and shows that this quantity can be evaluated analytically. Chapter 2 also addresses the estimation 
of error in the estimate, which is obtained using simulated catalogues. Other methods of 
obtaining this error, and also the bias in the maximum likelihood estimator, are outlined in the 
appendix. 


The results are presented in Chapter 3, which contains estimates of b value for the entire 
catalogue, temporal subsets of it, and for two subregions centered around Loppersum and Ten 
Boer. Chapter 4 contains the conclusions and recommendations for future work. 



SR. 15.BinnedParetoMLE 


-3- 


Restricted Draft Report 


2. Maximum Likelihood Estimator for Binned Seismic Data 

As mentioned in the introdnction, we take magnitude 1.5 to be the threshold of reliable 
observation. Since earthquake magnimdes are rounded to the nearest 1.5, this corresponds to a 
bin spanning a range of magnitudes 1.45 to 1.55. Multiplying Eq. (1.5) by tdn(lO), exponentiating 


and rearranging, we obtain 

M 0 — M c 10 dm ° (2.1) 

Setting m 0 =1.45, d = 1.5 and M c = 1.259 x 10 9 N-m in Eq. (2.1), we obtain 

M 0 =1.884x10" N-m (2.2) 

We shall work with the normalized moment 

x = M/M 0 (2.3) 

This means that the first bin for x starts at x = 1 and ends at CC, where 

cc = 10 ° u =10 015 =1.413 (2.4) 


The /èth bin runs from x — ot to ct. If the wth bin is the last non-empty bin, all the empty bins 
following it are included in the ;/th bin, so that it extends from x = Cl' ' to °o. The probability p k 


that an event falls within the ;/th bin is given by 

p k = F c (a k ~ 1 )-F c (a k )-k= 1,2 ... n- 1 (2.5a) 

p n =F c (a"- 1 ) (2.5b) 

where (gfEq. 1.6) 

F c (x) = x~ p (2.6) 

It follows that 

p k = cC P(k -" - a ~ pk ; k = 1, 2 1 (2.7a) 


(2.7b) 
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For a catalogue containing N events, the probability that bins 1, 2 to n contain q, r 2 to r n events is 
given by 


P(r 1 ,r 2 ...r n ) = N\fl{p k rk /rj} 

k=l 


(2.8) 


where ^ r A . = A (2.9) 

>t=i 


If we are given the r k and wish to make inferences abont the p k , the right-hand side of Eq. (2.8) 
can be expressed as a likelihood function Up x ,p 2 . ■ p „). That is 

L( Pl ,p 2 ...p n ) = N\Y[{p k rt /r k '} (2.10) 

k =1 


According to the maximum likelihood method, the best value of the parameter (5 controlling the 
p k that fits the observations (the values of is the one that maximizes L(p u p 2 ...p^). Assuming 
that L is a smooth function of j3, f5 maximizes L locally if 


dL 

dj3 


= 0 


( 2 . 11 ) 


and 



( 2 . 12 ) 


We now show that the value of f5 satisfying Eq. (2.11) can be determined analytically when the p k 
are given by Dividing Eq. (2.11) by E we may write it as 


d 

d(5 


ln(L) = 0 


(2.13) 


Substituting Eq. (2.10) into Eq. (2.13) we deduce that 


I 


k =1 


a 

¥ 


ln(p k ) = 0 


(2.14) 


Using Eq. (2.7a) we have 


ln (p k ) = -f5k ln(or) + ln(a^ — 1) ; h - 1 , 2 ... n -1 


(2.15a) 
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\n(p n ) = -jB(n-\)\n(a) 


(2.15b) 


so that — \ kr k —r n lln (a) + (N — r n ) ln {a)a^ /( a p - 1) = 0 


^=i 


(2.16) 


Eq. (2.16), and hence Eq. (2.11) is satisfied for 


1 -a p =- 


N - r 


2 >* 


This can be written as 


(2.17) 


cr fi = ^ 


k =1 _ 

n 

Z^ ~ r „ 


(2.18) 


Hence the value of j3 satisfying Eq. (2.18), which we henceforth denote by j3 MLE , is given by 


Pmle 


ln] Z kr k ~ r n f - ln ] tik- 1),, 

k =1 I 1 k =1 


ln(or) 


(2.19) 


Note that /? VIIË is infinite in the special case that r x — N and r 2 — r 3 — ... r N — 0. Hence we 

N 

exclnde this case when evaluating the statistics of /3 MlSE . Let '^ i f(r l ,r 2 ...r n ;e) denote the sum 

(r, ,-.r„} 

n 

of a function /(r,, r 2 ... r n ; £) over all non-negative integers r u r 2 ... r n snch that Z r k ~ N and let 

k=l 

N 1 

Z denote the same snm, but with the restriction that r 2 , r 3 ... r n are not all zero. We consider 
the restricted average of J3 KU I . given by 


M mle =<Pmle>'=N\ Z ‘X\(Pk rk !r k \) Z Pmle^-p") 


( 2 . 20 ) 


k= 1 

Similarly, the variance of /^ [LE is given by 


&’mle < Pmle ^ ^ Pmle 


(2.21) 
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where</^ 2 >'=iV! £ 'iW* /r * ! ) £ P^KX-p”) 

( r i .-'■»} *=1 


( 2 . 22 ) 


In this report we estimate // MLE and C^ [LE by performing Monte-Carlo simulations of catalogues 
with a specified /?. If u is a random number uniformly distribnted between 0 and 1, representing 
F(x), then Eq. (2.6) may be used to obtain At: 

jc = (1 -m)" 1//? (2.23) 


The value of 1 + INT{ln(.x)/ln(fl)}, where INT denotes integer part, then determines the 
nnmber of the bin in which x is located. 


As we shall see in the next chapter 1000 realizations of a catalogue is sufficiënt for confidence 
testing and error estimation. An alternative method of obtaining // MLE and Ö^ LE , based on 
transforming the summations on the right-hand sides of Eqs. (2.20) and (2.21) into integrations 
over one or two variables, is given in the appendix. These integrals may subsequendy be 
evaluated numerically or used to obtain large-IV asymptotics. 
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3. Results 

3.1 Maximum-Likelihood Estimate for the b-V alue of the Entire Groningen 
Catalogue and its Error Bounds 


The Groningen catalogue since l st May 1995 and up to the end of 2014 comprises 236 events of 
moment magnitude 1.5 or greater. The maximum observed magnitude was 3.6, observed for a 
single event. The seismic moments for the catalogue are therefore divided up into 22 bins, 
corresponding to moment magnitudes of 1.45-1.55, 1.55-1.65 and so forth, up to 3.45-3.55 and 
finally 3.55-°°. The number of events corresponding to each bin is tabulated in Table 3.1. 


Table 3.1. Distribution of Events in Seismic Moment Bins for the 236-Event Groningen 
Catalogue (l st May 1995-3l st December 2014, Rounded Moment Magnitude > 
1.5. 


Bin Number 

1 

2 

3 

4 

5 

6 

7 

8 

Moment Magnitude 

1.5 

1.6 

1.7 

1.8 

1.9 

2 

2.1 

2.2 

Number of Events 

54 

38 

26 

23 

18 

15 

7 

8 

Bin Number 

9 

10 

11 

12 

13 

14 

15 

16 

Moment Magnitude 

2.3 

2.4 

2.5 

2.6 

2.7 

2.8 

2.9 

3 

Number of Events 

8 

8 

8 

5 

4 

3 

1 

5 

Bin Number 

17 

18 

19 

20 

21 

22 



Moment Magnitude 

3.1 

3.2 

3.3 

3.4 

3.5 

3.6+ 



Number of Events 

0 

3 

0 

0 

1 

1 




The maximum likelihood estimate of the /?-value for this catalogue, according to Eq. (2.19) is 

Pmle= 0-644 (3.1) 

This corresponds to a maximum likelihood estimate of /Cvalue given by 

b MLE = 1.5x0.644 = 0.966 (3.2) 

This is so close to the accepted value b — 1 for induced seismicity that we use the value b - 1, 
corresponding to (3 — 2/3, to generate 1000 synthetic Monte-Carlo realizations of the 236-event 
Groningen catalogue, in order to obtain an error estimate for Eqs. (3.1) and (3.2). Figure 3.1 
shows the resulting distribution of /^i LE , which also indicates the simulation mean (0.664) and the 
actual catalogue value of /^, LE given in Eq. (3.1). The actual catalogue value is sufficiendy close to 
the simulation mean for the hypothesis that f3 — 2/3 to be accepted, which means that we can 
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use the Standard deviation of the simulations to provide an error bar on Eq. (3.1). This error bar 
is ± 0.043. Hence we may write 

Pmle - 0-644 ± 0.043 (3.3) 

The corresponding result for £ MLE is 

b MLE =0.966 ±0.065 (3.4) 



F/gure 3.1. Distribution Function of /? MLE obtainedfrorn 1000 Reali^ations of a 236-Event Catalogue with 

2/ 3. The vertical solid line shows the value of the simulation mean while the dashed line is 
the actual Groningen Catalogue /? MLE . 


The simulation is of value in providing further information about the statistics of the event 
catalogues. For example, a lower [3 value results in a greater probability of higher magnitude 
events, so we can ask the question about how well the event, or events, of maximum magnitude 
in a simulated catalogue correlate with the /4 ILE f° r that catalogue. The results are shown in 
Figure 3.2. The highest magnitude event, a 6.9, occurs for a /^ ILE of 0.629, only a litde lower than 
the simulation mean, while of the three catalogues for which the highest magnitude event is only 
magnitude 3, two of them have high values of /? MLE , 0.747 and 0.750, while the remaining one has 
a /? MLE of only 0.612. The complementary distribution function for the maximum magnitude 
event, or events, in a catalogue, is shown in Figure 3.3. We take note of the fact that 79% of the 
simulated catalogues contain an event of maximum magnitude greater than 3.6, the observed 
value for the actual catalogue. The shape of the curve in Figure 3.3 may be readily understood by 
considering the probability P N (x) that at least one event out of N has a moment magnitude 
greater than x. This is given by 
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P,W = l-{l-i ? cW} ï (3-5) 

where F c (x) is given by Eq. (2.6). The probability that at least one event has moment magnitude 
greater than m, P N (m), is then given by 

P N (m) = P N (x) (3.6) 

where x = 70 ( '"“ L45) (3.7) 

Substituting Eq. (3.6) into Eq. (3.5), and applying Eqs. (3.7) and (2.6), we obtain 

P N (m) = 1 - {1 - io- /J(m - | - 45, } A ' (3.8) 


Eq. (3.8) is also included in Figure 3.3 and is in close agreement. 


Maximum Magnitude versus LE for 1000 
Realizations of the 
236-Event Groningen Catalogue 
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Figure 3.2. Maximum magnitude event(s) in simulated catalogue versus /i MLE . 
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Figure 3.3. 


3.2 


Complementary Distribution Function of Maximum Magnitude 
Events from 1000 Realizations of a Simulated Groningen Catalogue 



Maximum Magnitude Event(s) m in Catalogue Realization 


Complementary Distribution Function of Maximum Magnitude Events from 1000 
Realizations of a 236-Event Groningen Catalogue. The Red Curve is the Analytical 
Expression in Eq. (3.8) 

Maximum Likelihood and Error Estimates for the />Value of Temporal 
Subsets of the Groningen Catalogue 


In this section the Groningen catalogue is divided into four sub-catalogues, each containing 59 
events, according to temporal order. 
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Table 3.2. Distribution of Events in Seismic Moment Bins for the 236-Event Groningen 
Catalogue, Split Up by Temporal Sub-catalogue 


Bin Number 

1 

2 

3 

4 

5 

6 

7 

8 

Moment Magnitude 

1.5 

1.6 

1.7 

1.8 

1.9 

2 

2.1 

2.2 

# Events T* Subcatalogue 

11 

16 

6 

6 

3 

3 

2 

2 

# Events 2 nd Subcatalogue 

13 

6 

7 

5 

5 

3 

2 

5 

# Events 3 rd Subcatalogue 

16 

10 

8 

7 

5 

3 

0 

1 

# Events 4 th Subcatalogue 

14 

6 

5 

5 

5 

6 

3 

0 

Bin Number 

9 

10 

11 

12 

13 

14 

15 

16 

Moment Magnitude 

2.3 

2.4 

2.5 

2.6 

2.7 

2.8 

2.9 

3 

# Events T* Subcatalogue 

1 

2 

2 

1 

2 

0 

0 

2 

# Events 2 nd Subcatalogue 

4 

1 

2 

3 

0 

0 

0 

1 

# Events 3 rd Subcatalogue 

2 

2 

4 

0 

0 

0 

0 

0 

# Events 4 th Subcatalogue 

1 

3 

0 

1 

2 

3 

1 

2 

Bin Number 

17 

18 

19 

20 

21 

22 



Moment Magnitude 

3.1 

3.2 

3.3 

3.4 

3.5 

3.6+ 



# Events T‘ Subcatalogue 

0 

0 

0 

0 

0 

0 



# Events 2 nd Subcatalogue 

0 

1 

0 

0 

1 

0 



# Events 3 rd Subcatalogue 

0 

1 

0 

0 

0 

0 



# Events 4 th Subcatalogue 

0 

1 

0 

0 

0 

1 




The maximum likelihood estimates of the j3-v alue for each sub-catalogue, according to Eq. (2.18) 
are (with corresponding Avalue maximum likelihood in brackets): 


l st sub-catalogue 

Pmle 

= 0.695 (b MLE 

= 1.043) 

(3.5a) 

2 nd sub-catalogue 

Pmle 

= 0.577 (V, : 

= 0.865) 

(3.5b) 

3 rd sub-catalogue 

Pmle 

= 0.809 (b MLE 

= 1.213) 

(3.5c) 

4 th sub-catalogue 

Pmle 

= 0.528 (b MLE 

= 0.782) 

(3.5d) 


A simulation of 1000 realizations of a 59-event catalogue with ƒ? — 2/3 yields the result 


( 3 mle =0.665 ±0.086 


(3.6) 
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Inspection of Eqs. (3.5) reveals that all the sub-catalogues with the exception of the fïrst lie 
outside the error bars in Eq. (3.6). This is represented graphically in Figure 3.4, where the values 
of Aile given in Eqs. (3.5) are compared with the distribution of /? MLE obtained from the 
simulation. 


Distribution Function of Maximum Likelihood Estimate 
from 1000 Realizations of a 59-Event Catalogue with f3— 2/3 
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Figure 3.4. Distribution Function of Maximum Fikelihood Estimator Obtainedfrom 1000 FeaUgations 
of a 59-Event Catalogue with f3 — 2f3. The Simulation Mean, Together with the Four 
Values of f3 MLE Given in Eqs. 93.45), are Displayed as Vertical Lines on the Vlot. 

We see from Figure 3.4 that the hypothesis that the f3 value for the fourth sub-catalogue is less 
than 2/3 has a confidence level greater than 95%, and that we can be almost as confident that the 
f3 value for the third sub-catalogue is greater than 2/3. Note that the probability that at least two 
of the sub-catalogues have /? MLE values in the bottom or top 5% of the distribution is 1 — 0.9 4 — 4 
x 0.1 x 0.9' = 0.0523. We now determine whether there is any correlation with rate effect by 
examining the time intervals corresponding to the four sub-catalogues. The dates and times of 
the fïrst and last events in each sub-catalogue are listed in Table 3.3. The maximum likelihood 
estimate for each sub-catalogue is plotted against the catalogue duration in Figure 3.5. It is hard 
to draw any conclusions regarding rate dependency from Figure 3.4, and we further subdivide the 
catalogue into la, lb, 2a, 2b etc. where the a- half of each of the four sub-catalogues consists of 29 
events and the £-half 30 events. The number of events by magnitude for each of these sub- 
catalogues is listed in Table 3.4. 


Table 3.3 Dates and Times of First and Tast Events for Each of the Four 59-Event Temporal Sub- 
Catalogues of the 236 Event Groningen Catalogue 



Date 

Time 

Date 

Time 

Elapsed Hrs. 

l st Sub-Catalogue 

15/05/1995 

08:24:00 

25/02/2005 

20:38:24 

3586.24 

2 nd Sub-Catalogue 

10/03/2005 

18:43:12 

08/05/2009 

22:48:00 

1524.08 

3 rd Sub-Catalogue 

05/07/2009 

08:52:48 

19/02/2012 

02:52:48 

953.00 

4 th Sub-Catalogue 

31/03/2012 

09:21:36 

30/12/2014 

Not Recorded 

1006.64* 


Time of last event for fourth sub-catalogue is assumed to be 12:00 pm 
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Maximum Likelihood Estimator versus Duration for 59-Event 
Temporal Sub-Catalogues of the Groningen Catalogue 
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Figure 3.5. Maximum Uke/ihood FLstimators of b-Value forFour 59-Event Sub-Catalogues of the 

Groningen Catalogue versus Sub-Catalogue Duration. 


The maximum likelihood estimate has been computed, according to Eq. (2.19), for each of the 
sub-catalogues in Table 3.4, and these have been plotted against the duration of the cub- 
catalogue in Figure 3.6. For these sub-catalogues containing almost equal number of events (29 
or 30), it is now clear that there is little correlation between the maximum likelihood estimate of 
the parameter f3 and the duration of the sub-catalogue. The error bars (corresponding to ± one 
Standard deviation) shown in Figure 3.6 are obtained from simulations of 1000 realizations of 
each sub-catalogue with the f3 for the underlying process set equal to /? MLE . 
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Table 3.4. Distribution of Events in Seismic Moment Bins for the 236-Event Groningen 
Catalogue (l st May 1995-3l st December 2014, Ronnded Moment Magnitude > 
1.5), Split Up into Eight 29 or 30 Event Temporal Sub-catalogues 


Bin Number 

1 

2 

3 

4 

5 

6 

7 

8 

Moment Magnitude 

1.5 

1.6 

1.7 

1.8 

1.9 

2 

2.1 

2.2 

# Events 1 a Subcatalogue 

4 

9 

2 

5 

1 

3 

1 

0 

# Events 1 b Subcatalogue 

7 

7 

4 

1 

2 

0 

1 

2 

# Events 2a Subcatalogue 

7 

5 

2 

3 

3 

1 

0 

2 

# Events 2b Subcatalogue 

6 

1 

5 

2 

2 

2 

2 

3 

# Events 3a Subcatalogue 

9 

5 

4 

3 

0 

2 

0 

0 

# Events 3b Subcatalogue 

7 

5 

4 

4 

5 

1 

0 

1 

# Events 4 a Subcatalogue 

7 

1 

3 

4 

2 

4 

0 

0 

# Events 4 b Subcatalogue 

7 

5 

2 

1 

3 

2 

3 

0 

Bin Number 

9 

10 

11 

12 

13 

14 

15 

16 

Moment Magnitude 

2.3 

2.4 

2.5 

2.6 

2.7 

2.8 

2.9 

3 

# Events 1 a Subcatalogue 

0 

1 

2 

1 

0 

0 

0 

0 

# Events 1 b Subcatalogue 

1 

1 

0 

0 

2 

0 

0 

2 

# Events 2a Subcatalogue 

3 

1 

1 

0 

0 

0 

0 

0 

# Events 2b Subcatalogue 

1 

0 

1 

3 

0 

0 

0 

1 

# Events 3a Subcatalogue 

2 

2 

2 

0 

0 

0 

0 

0 

# Events 3b Subcatalogue 

0 

0 

2 

0 

0 

0 

0 

0 

# Events 4 a Subcatalogue 

0 

3 

0 

0 

2 

0 

0 

1 

# Events 4 b Subcatalogue 

1 

0 

0 

1 

0 

3 

1 

1 

Bin Number 

17 

18 

19 

20 

21 

22 



Moment Magnitude 

3.1 

3.2 

3.3 

3.4 

3.5 

3.6+ 



# Events 1 a Subcatalogue 

0 

0 

0 

0 

0 

0 



# Events 1 b Subcatalogue 

0 

0 

0 

0 

0 

0 



# Events 2a Subcatalogue 

0 

0 

0 

0 

1 

0 



# Events 2 /j Subcatalogue 

0 

1 

0 

0 

0 

0 



# Events 3a Subcatalogue 

0 

0 

0 

0 

0 

0 



# Events 3b Subcatalogue 

0 

1 

0 

0 

0 

0 



# Events 4a Subcatalogue 

0 

1 

0 

0 

0 

1 



# Events 4ó Subcatalogue 

0 

0 

0 

0 

0 

0 
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Maximum Likelihood Estimator versus Duration for 29&30-Event 
Temporal Sub-Catalogues of the Groningen Catalogue 
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F/gure 3.6. Maximum Fikelihood Fstimators of b-J^alueforEight 29 or 30-Fvent Sub-Catalogues of the 

Groningen Catalogue versus Sub-Catalogue Duration. 

Further ways of temporally subdividing the Groningen catalogue into equal or nearly equal 
subcatalogues are listed in Table 3.5, and the resulting maximum likelihood estimators of (3 are 
shown in Figure 3.7. These results are compared with the hypothesis that the underlying value of 
/?for each sub-catalogue is in fact equal to 2/3 (corresponding to b — 1) in Figure 3.8. The 
horizontal lines represent the limits of 95% confidence intervals for this hypothesis, obtained 
from simulations containing 5000 realizations. 


Table 3.5. Subdivision of the 256-Event Groningen Catalogue into up to Eight Sub-Catalogues 


Number of Sub-Catalogues 

Subdivision 

1 

1 x 236 

2 

2 x 118 

3 

79, 78, 79 

4 

4 x 59 

5 

4 x 47, 48 

6 

40, 4 x 39, 40 

7 

33, 5 x 34, 33 

8 

29, 30, 29, 30, 29, 30, 29, 30 


We see that the computed values of /^ ILE are consistent with the hypothesis that /? = 2/3 at the 
95% level for all but four of the sub-catalogues, and these four are on the lower bound of the 
confidence interval rather than being significantly below it. 
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Maximum Likelihood Estimates of /^-Parameter for Sub- 
Catalogues of the Groningen Catalogue 
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Figure 3. 7. Maximum Likelihood Estimators of (3for Sub-Catalogues of the 256-Event Groningen Catalogue 
Selected as Shown in Table 3.5. 
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Figure 3.8. 


As Figure 3 .7 but als o showing the 95% Confidence Intervals for the Hypothesis that 13 — 2/3 for 
Each Sub-Catalogue. Note that the Smaller the Sub-Catalogue the Wider the Confidence Interval. 
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3.3 Maximum Likelihood and Error Estimates for the b-Value of Events in the 
Loppersum and Ten Boer Areas 

We now focus on two specific spatial sub-regions of the Groningen field centered around 
Loppersum (containing the region of historically the most intense seismicity, and also the event 
of maximum moment magnitude) and Ten Boer, the most active region outside Loppersum, and 
also closer to the city of Groningen. These regions have centres (244, 598) and (250, 591) 
respectively, in terms of Easting and Northing coordinates expressed in kilometers. The radius of 
each region is 5 kilometers. These regions correspond to the ones chosen in [5], and the estimates 
of ^-value obtained here can be compared with those obtained in [5]. The number of events of 
magnitude 1.5 or greater in each region during the period l st May 1995 to 31 st December 2014 is 
82 for Loppersum and 60 for Ten Boer. The events are tabulated by moment magnitude in Table 

3.4. The spatial location of these events, as well as an indication of magnitude, is shown in Figure 

3.5. 
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Table 3.6. Distribution of Events in Seismic Moment Bins for the Loppersum and Ten Boet 


regional sub-catalogues (l st May 1995-3T' December 2014, Ronnded Moment 
Magnitude > 1.5). 


Bin Number 

1 

2 

3 

4 

5 

6 

7 

8 

Moment Magnitude 

1.5 

1.6 

1.7 

1.8 

1.9 

2 

2.1 

2.2 

Loppersum Region 

14 

8 

7 

6 

7 

8 

2 

6 

Ten Boer Region 

15 

10 

11 

4 

4 

4 

0 

1 

Bin Number 

9 

10 

11 

12 

13 

14 

15 

16 

Moment Magnitude 

2.3 

2.4 

2.5 

2.6 

2.7 

2.8 

2.9 

3 

Loppersum Region 

4 

2 

3 

2 

4 

1 

1 

3 

Ten Boer Region 

1 

4 

2 

1 

0 

0 

0 

2 

Bin Number 

17 

18 

19 

20 

21 

22 



Moment Magnitude 

3.1 

3.2 

3.3 

3.4 

3.5 

3.6+ 



Loppersum Region 

0 

2 

0 

0 

1 

1 



Ten Boer Region 

0 

1 

0 

0 

0 

0 




Note that the events are color-coded according to whether they occur in the first half (by time of 
event) or second half of the sub-catalogue. These two halves are considered separately in the next 
section. 

Applying Eq. (2.18) to the data in Table 3.5, we obtain the following values of Ale for the 
Loppersum and Ten Boer sub-catalogues (the corresponding values of £ MLE are given in 
brackets): 


Loppersum 

Ale = 0.471 (£ mle = 0.707) 

(3.7) 

Ten Boer 

Ale = 0.709 (£ mle = 1.064) 

(3.8) 


1000 realizations of the 82-event Loppersum sub-catalogue and the 60-event Ten Boer 
sub-catalogue with input value /? = 2/3 yielded Au; = 0.668 ± 0.074 for Loppersum and 0.669 ± 
0.085 for Ten Boer. The actual value of Ale is far outside the error bars of the simulation for 
Loppersum and well within it for Ten Boer. Further illustration of this is provided by Figure 3.6 
which shows the distribution of Ale for the two simulations together with the actual values in 
Eqs. (3.7) and (3.8). The shapes of the distributions are compared with Gaussian distributions 
with the same mean and variance for comparison. Figure 3.6 (a), in particular, tells us that we can 
be 99.9% confident that the underlying process producing induced seismicity in the Loppersum 
region has a value smaller than 2/3. Note that the estimate of /?obtained by the present 
method for the Loppersum sub-region is smaller than that obtained in [5] using a straight line fit 
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of log frequency versus seismic moment. The value obtained there was f5 = 0.57, corresponding 
to b — 0.85. This method for estimating ^-value has been used in the past, but as pointed out in 
[6], is not justified on statistical grounds. 


Spatial, Temporal and Magnitude Distribution of Events in the 



Distance East of Loppersum Centre (km) 


Spacial, Temperal and Magnitude Distribution of Events 
in the Ten Boer Sub-Catalogue 


• m = 3.0-3.2, lst Half 
O m = 2.5-2.9, lst Half 
O m = 2.0-2.4, lst Half 
o m = 1.5-1.9, lst Half 

• m = 3.0-3.2, 2nd Half 
O m = 2.5-2.9, 2nd Half 
O m = 2.0-2.4, 2nd Half 
o m = 1.5-1.9, 2nd Half 

Bounding Circle 

-6 -4 -2 0 2 4 6 



Distance East of Ten Boer Centre (km) 


F/gure 3.9. Spatial distribution of events in (a) the Toppersum and (b) Ten Boer regions. Magnitude classes 
are also indicated. The terms V' Half and 2 tul Half refe r to the first and second halves of the 
catalogues, arranged temporalij 
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Figure 3.10. Distribution of maximum likelihood estimate /? M , , for 1000 realfations of 
simulated (a) Foppersum and (b) Ten Boer sub-catalogues mth ( 3 — 2/3. The 
shapes of the distributions are compared mth Gaussian distributions having the 
same mean and variance and the values of /? MLE obtained for the actual sub- 
catalogues are also shonm. 
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3.4 Maximum Likelihood and Error Estimates for the b-V alue of Temporal 

Subsets of Events in the Loppersum and Ten Boer Areas 


Because the number of events in the Loppersum and Ten Boer sub-regions is relatively small, 
just two temporal subsets, corresponding to the first and second halves of the sub-catalogue 
corresponding to each sub-region, are used. The distribution of the events into the magnitude 
bins for each temporal subset is shown in Table 3.6. 

Table 3.7. Distribution of Events in Seismic Moment Bins for the first and second halves of 


the Loppersum and Ten Boer subcatalogues (l st May 1995-3E' December 2014, 
Rounded Moment Magnitude > 1.5). 


Bin Number 

1 

2 

3 

4 

5 

6 

7 

8 

Moment Magnitude 

1.5 

1.6 

1.7 

1.8 

1.9 

2 

2.1 

2.2 

Loppersum Region (T‘ Half) 

7 

7 

1 

2 

3 

5 

1 

2 

Loppersum Region (2 nd Half) 

7 

1 

6 

4 

4 

3 

1 

4 

Ten Boer Region (T‘ Half) 

7 

7 

6 

1 

2 

1 

0 

1 

Ten Boer Region (2 nd Half) 

8 

3 

5 

3 

2 

3 

0 

0 

Bin Number 

9 

10 

11 

12 

13 

14 

15 

16 

Moment Magnitude 

2.3 

2.4 

2.5 

2.6 

2.7 

2.8 

2.9 

3 

Loppersum Region (T‘ Half) 

4 

2 

2 

1 

2 

0 

0 

1 

Loppersum Region (2 nd Half) 

0 

0 

1 

1 

2 

1 

1 

2 

Ten Boer Region (l st Half) 

0 

2 

1 

1 

0 

0 

0 

1 

Ten Boer Region (2 nd Half) 

1 

2 

1 

0 

0 

0 

0 

1 

Bin Number 

17 

18 

19 

20 

21 

22 



Moment Magnitude 

3.1 

3.2 

3.3 

3.4 

3.5 

3.6+ 



Loppersum Region (T‘ Half) 

0 

0 

0 

0 

1 

0 



Loppersum Region (2 nd Half) 

0 

2 

0 

0 

0 

1 



Ten Boer Region (T‘ Half) 

0 

1 

0 

0 

0 

0 



Ten Boer Region (2 nd Half) 

0 

0 

0 

0 

0 

0 




The values of /? vll l , for each temporal sub-catalogue, computed using Eq. (2.19), are listed in Eqs. 
(3.9a)-(3.9d), with the corresponding values of b MlSE in brackets. 


Loppersum, l st half 

^=0.500(^=0.751) 

(3.9a) 

Loppersum, 2 nd half 

P MLE = 0-435 (b MLE - 0.652) 

(3.9b) 



SR. 15.BinnedParetoMLE 


-22- 


Restricted Draft Report 


Ten Boer, l st half 

f3 MLE =0.764 (b MLE =1.146) 

(3.9c) 

Ten Boer, 2 nd half 

Pmle =0-661 (b MLE = 0.992) 

(3.9d) 


In order to test whether there is a significant temporal shift in the values of /? MLE for the 
Loppersnm and Ten Boer sub-regions, simulations of 41 and 30 event catalogues were 
performed, these being the nnmber of events in the Loppersum and Ten Boer temporal sub- 
catalogues, respectively. The values of /3 used as input were the maximum likelihood estimates of 
the entire catalogue for each subregion, namely the values given in Eqs. (3.7) and (3.8) 
respectively. The results are shown in Figure 3.7. 


Distribution Function of Maximu Likelihood Estimate from 
1000 Realizations of a 41-Event Catalogue with f3 = 0.471 



0 0.2 0.4 0.6 0.8 

Maximum Likelihood Estimate 


Distribution Function of Maximu Likelihood Estimate from 
1000 Realizations of a 30-Event Catalogue with p— 0.709 



Maximum Likelihood Estimate /? V1 T1 , 


Figure 3.9. Distribution of maximum likelihood estimate f3 WL ^for 1000 realizations of simulated (a) 

Toppersum and (b) Ten Boer sub-catalogues with f3- 2/3. The shapes of the distributions are 
compared with Gaussian distributions having the same mean and variance and the values of 
Aile obtained for the actual sub-catalogues are also shown. 
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We can conclude from Figure 3.6 that the differences in the /? MLE observed between the first and 
second halves of the catalogues for Loppersum and Ten Boet are not significant enough to 
indicate a change in the underlying process. 
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4. Conclusions and Recommendations 

This report contains a prescription for obtaining a maximum likelihood estimate for the Avalue 
in the Gutenberg-Richter relationship for the frequency versus moment magnitude for seismic 
events. This relationship has been applied for both naturally occurring and induced seismicity. A 
procedure for obtaining the error estimate in the resulting values of b is also outlined. 


The prescription and procedure is then applied to the catalogue of seismic events triggered by gas 
production in the Groningen field. Events occurring between l st May 1995 and 31 st December 
2014, and of magnitude 1.5 and greater are used to estimate values of b for the entire catalogue 
and various subsets of it. The value of b obtained for the entire catalogue is consistent with the 
commonly accepted value of unity for induced seismicity. 


A key question for hazard assessment is whether the character of the seismicity is changing with 
time, or is dependent on the number of events per unit time. This is because the frequency of 
rare events of high magnitude, potentially capable of causing structural damage, is dependent on 
the Avalue. Accordingly, the Groningen catalogue is then applied to temporal subsets to examine 
any evidence for a varying Avalue. While a subdivision of the catalogue into four subsets of equal 
numbers of events suggests more variability in the maximum likelihood estimator for the subsets 
than would be expected if the Avalue of the underlying process were constant, a closer 
examination, involving further subdivisions of the catalogue, does not show any systematic 
dependence of b value on event rate. 

A recent report examined induced seismicity in various spatial sub-regions of the Groningen field 
and beyond it and related the data to a stochastic model of fault-slippage induced seismicity. Two 
of these sub-regions, around Loppersum, where the event of largest magnitude to date has 
occurred, and around Ten Boer, which is closer to the city of Groningen, were selected for 
analysis by the methods introduced in this report. All events of magnitude greater or equal to 1.5 
and occurring between l st May 1995 and 31 st December 2014 were considered. The results 
showed that the Avalue around Loppersum is less than unity to a high degree of significance, 
while the Ten Boer data is consistent with a Avalue of unity. No statistically significant temporal 
evolution of Avalue for these spatial sub-regions was found. Note that no attempt to obtain a 
spatial correlation dimension for events, such as was done for seismicity in the Lacq gas field [7], 
was made here. 

Though it resembles the one introduced in [7], the method presented here for determining A 
value is different from any proposed in the literature and it should be applied to literature 
datasets in order to compare with these other methods. Furthermore, the Avalue estimator is 
largely dependent on a single parameter, the summation of the number of events occurring in a 
magnitude moment bin times the bin number. Examination of the magnitude of the likelihood 
function itself should be performed in order to identify data anomalies or inconsistencies. 



SR. 15.BinnedParetoMLE 


- 25 - 


Restricted Draft Report 



SR. 15.BinnedParetoMLE 


-26- 


Restricted Draft Report 


5. Nomenclature 

Roman Characters 


Symbol 

Definition 

Unit 

Introduced 

b 

Avalue in Gutenberg-Richter law 

- 

Page II 

^MLE 

Maximum likelihood estimator of b 

- 

Page II 

d 

Constant relating moment and moment magnitude 

- 

Eq. (1.1) 

F(M) 

Distribution function for seismic moment 

- 

Eq. (1.6) 

Fc 

Complementary distribution function 

- 

Eq. (2.6) 

fir i-0 

Arbitrary function of bin event numbers 

- 

Eq. (2.20) 

Upi-pj 

Likelihood function 

- 

Eq. (2.10) 

M 

Seismic moment 

N-m 

Eq. (1.1) 

Mo 

Threshold seismic moment 

N-m 

Eq. (2.1) 

M c 

Reference seismic moment 

N-m 

Eq. (1.1) 

m 

Moment magnitude 

- 

II 

m 

Threshold moment magnitude 

- 

II 

N 

Number of events in a catalogue 

- 

Eq. (2.8) 

n 

Largest bin number containing an event 

- 

Eq. (2.5) 

P(m) 

Probability that an event has moment magnitude 

- 

II 

-Pn(/w) 

Probability that the largest event in an N-event catalogue has 
moment magnitude > m 

- 

Eq. (3.6) 

P(n..r n ) 

Probability that rk events fall into the kth bin for all bins 1 to n 

- 

Eq. (2.8) 

pk 

Probability that an event falls in dre /èth bin 

- 

Eq. (2.5) 

n. 

Number of events in the kth bin 

- 

Eq. (2.8) 

So 

Component of maximum likelihood estimator 

- 

Eq. (A.2) 

Si 

Component of maximum likelihood estimator 

- 

Eq. (A.2) 

Si(è) 

Generalization of Ti 

- 

Eq. (A.9) 

u 

Uniformly distributed random variable over [0,1] 

- 

Eq. (2.23) 

X 

Normalized seismic moment ( M/Mo ) 

- 

Eq. (2.3) 
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Greek Characters 


Symbol 

Definition 

Unit 

Introduced 

CC 

Ratio of moments for a moment magnitude increment of 0.1 

- 

Eq. (2.4) 

fi 

Exponent in Pareto distribution of nromalized moment 

- 

Eq. (1.7) 


Maximum likelihood estimator of 

- 

Eq. (2.19) 

£ 

Small quantity set to zero at end of calculation 

- 

Eq. (A.9) 

/4ile 

Mean of distribution of /Ui.e 

- 

Eq. (2.20) 


Variance of distribution of /?hle 

- 

Eq. (2.21) 
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Appendix 1 — Conversion of the Summations in Eqs. (2.20) & (2.22) into 
Integrations 

In this appendix a methodology for converting the summations on the right-hand sides of Eqs. 
(2.20) and (2.22) into integrations over a single variable is outlined. Considering Eq. (2.19) first of 
all, we write 


Pmle ={^o - Aj/lnCör) 


(A.1) 


where S 0 = ln] ^ kr k - r B i; S l = ln] ^ (k - 1> 


(A.2) 


The conversion is accomplished by using the following integral representation of the logarithm: 


ln(z) = J duu 1 {exp(-n) - exp(-zw)} 


(A.3) 


Hence S 0 = j du\ exp (-u) - expj - ^ kr k 


\k= 1 


>u >/u 


(A.4) 


The unrestricted summation over T 0 may now be evaluated explicitiy on reversing the order of 
summation and integration to yield 


< S 0 >=jduu ' \exp(-u) -(^e uk p k +e l) p„) N [ (A.5) 

o l *=i J 

The restricted sum may readily be related to the unrestricted one by subtracting the term ln(N)/> 1 N 
and dividing by 1 — p t N . The result is 


<S 0 > , = {<S 0 >-ln(N)p l N }/(l-p l N ) 


(A.6) 


Integrating Eq. (A.5) by parts, we obtain 


n —1 


< S 0 >— 1 du ln(n)] exp(-n) -l- e uk p k +e u(n l) p n ) r 

o l du k=1 

So that we finally have 


(A.7) 


<s 0 >'=*-' HN)P ‘ 


1 N j du ln(n) jexp(-n) + (^e uk p k +<? l) p n ) r 


(A.B) 



SR. 15.BinnedParetoMLE 


-30- 


Restricted Draft Report 


We proceed in a slightly different way for S\ due to the fact that <i’ 1 > is infinite. We introducé 
the quantity 


+ X (*-!)/ 


where fis a positive quantity that we set to zero at the end of the calculation. We now have 


< S,(f) >' = {< S^e) >-ln(e) Pl N }/(l- Pl N ) 


(A.10) 


Using Eq. (A.3), we have 


S 1 (£) = j du\ exp(-w) - expj - £ + ^\(k -1 )r k \u l/ 


(A.11) 


Interchanging the order of integration and summation in the unrestricted sum, and evaluating the 
sum yields 


< 5; (f) >= j dun' 1 < exp(-w) - exp(-a/) ^ Pn exp{-(& -1 )u } 


Since ln(f) = J duu 1 {exp(-i«) - exp(-a*)} 


(A.12) 


(A.13) 


we may wnte 


< S^e) >'= jduu 1 1 exp(-w)-exp(-a/H ^ Pk exp{-(A: -1 )u) - Pl N >/(l — Pl N ) 


(A.14) 


Integration by parts yields the result 

“ f d (n \ N 

< 5, {£) >' = j du\n(u)< exp(-u) + exp(-at) — ^ Pk expj-i/t -1);/} /( 1- Pl N ) 

o du\ k= i ) 


f|t/idn(M)exp(-a<)j ^ Pk exp{-(fc - \)u\ - Pl N 1/(1 — Pl N ) 


(A.15) 


As £—> 0 the first term on the right-hand side of Eq. (A.15) tends to a constant and the second 
term tends to zero as £ Hence, 


< Sj >=< SjCO) >' = |du (1- Pl N )exp(-M) + — exp{-(£ -l)w} l (A.16) 
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For numerical integration purposes, or to develop asymptotic expansions for large N, it is useful 
to transform to the new variables v x and v 2 in Eqs. (94) and (95) respectively as follows: 


exp(-Vj) = J^e~ uk p k +e~ u(n ~ l) p n 


(A.17) 


and exp(-v 2 ) = ^ e “ <k 1+£) p k 


(A.18) 


Raising Eqs. (A.17) and (A.18) to the power N and differentiating, we obtain 


d\\N exp(-\\N) =—— (Ye " k p t +e " ( " l) p„) N du 
du 


(A.19) 


and t/v 2 Aexp(-v 2 A) = ——(Y« “ ( * Ue) p k ) N du 

du fzf 


(A.20) 


On substituting Eqs. (98) and (99) back into Eqs. (94) and (95), we get 


< 


S 0 > = -— ln(jV) / 1 |j du\n(u)exp(-u) - ivj dv t ln(M)exp(-Av 1 )| 


(A.21) 


and < 5j > / = Jt/Mln(M)exp(-M)-A|c?v 2 ln(M)exp(-Av 2 )/(l- 1\) (A.22) 

o o 

Developing expansions for u in powers of v x in Eq. (A.21) and v 2 in Eq. (A.22) then leads to large 
N asymptotic expansions for <i’ 0 > / and <i’ 1 > / . 


We shall not deal in detail here with the representation of </? MLE 2>/ as an integral, only 
mentioning that we make nse of the fact that (cfEq. A.3) 

ln(z)ln(z / ) = jduu~ l Jt/vv _I {exp(—M)-exp(-zw)}{exp(-v)-exp(-z / v)} 
o o 


(A.23) 
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